#define _NOFMC
#define ARRTEST_
subroutine MCrun(nbatch, kinit, lact)
    use util_mc
    use define
	use allocs
    use statistics
!    use util_cmfd,      only : cmfd
!    use HMGZ
    use util_mc
    use RNG
!    use FMC
!    use THFDB,          only :fdb_th 
    use cmfd2g,         only :cmfd2gacc
    use lcmfd2g,        only :lcmfd2gacc
#ifdef MPIF
	use MPIComm_mod
#endif
	implicit none
    include "mcbase.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "accm.FI"
    include "xsec.FI"
    include "files.FI"
    include "FMC.FI"
    include "hmgz.FI"
    include "statistics.FI"
    include "GroupConstant.FI"
    include "mpi_arrays.FI"
!$    include "omp_lib.h"
!$    include "mtbuf.FI"
    integer, intent(in) :: nbatch
    real(8), intent(in) :: kinit
    logical, intent(in) :: lact
    integer :: base, g, g2, i, mi, j, k, l, r, t, jc, ic, gp, batch, iter
    integer, save :: nsource
    integer, save :: totbat=1
    real(8) :: chk, test, err, shnEtrp, shnEacc, shnECMFD
    real(8), save :: wtsum, k_mc, k_cmfd, k_cmfdp, k_nn, ke
    integer, save :: nnt
    integer :: batwp, ntally, nc_super
    logical :: lmwnd, loncmfd, lonfdb=.false.
    real(8) :: ts(2)
    integer :: tid
    integer :: nhtcnt(16)
    integer(8) :: lseed
    real(8) :: lkomp, nrmj, nrmp, volflux , phir, vcmr, scmrx, scmry
    integer :: nqompt(3) , tpin, nr
    real(8) :: tmc, tmcacc
    logical, save :: lfirstcmfd = .true.
!    integer :: addntrn
    real(8) :: rphiv
    real(8) :: tpsi, nrmpsi, psisumacc(1:32), psisumcyc(1:32)
    real(8) :: err_cmfd, err_cmfd_max, err_asmp, err_p2bf
    integer, parameter :: fpsi = 1010
    character(256) :: sfile
    real(8) :: leak, nrmv, nfis, nsct, lhs, rhs, diff, upper, lower, errblc
    real(8) :: lsupacc(32), lsupcyc(32), lslw(32), fpsidiff, perrcyc, perracc, fpsierr,rnormpsiacc, rnormpsicyc
    real(8), dimension(nfmx, nfmy) :: acpsi, acpsip
    logical, save :: once = .true.
    integer :: qi, qj, ai, aj, npxm, npxp, npym, npyp, fi, fj
    real(8) :: tilt_c, tilt_a, fss, fssa, qtp(2,2), qtpa(2,2), tiltx, tilty
	real(8) :: gkrcv, maxpinerr(32), pinerr
	integer :: nqtot, ia, ja, stx, sty
	integer :: nfsrnxt
	
	real(8) :: shnREC, rmserr_rec
	integer :: nwt
	integer, save :: mccount = 0
	real(8) :: s1, s2
	
	
	fpsidiff =0
	nfsrc=0
	nfsrcf=0

	ppsif=>psif
	ptot=> GCAtot
	pnuf=> GCAnuf
	pphiv=>GCAphiv
	psct=>GCAsct
	pjc=>hjac
#ifdef MPIF	
!	call dmalloc(psifrcv, nfmx, nfmy)
	ppsif=>psifrcv
	ptot=>totrcv
	pnuf=>nufrcv
	pphiv=>phivrcv
	psct=>sctrcv
	pjc=>jcrcv
#endif
    acpsi=0;acpsip=0;
    lsupacc=0;lsupcyc=0;lslw=0
    err_asmp=0; err_p2bf=0;
    
    if (RSTDHAT .and. lact) then
        GCAtot =0
        GCAnuf =0
        GCAphiv=0
        GCAsct =0
        hjac   =0 
!	    ptot =0
!	    pnuf =0
!	    pphiv=0
!	    psct =0
!	    pjc  =0
        logstring = "Reset Dhat."
        call log_(logstring)
    endif
    
!10  format     (i6,1x, i10,1x,  f7.5,1x, f8.3,1x, f8.3,1x, i3,1x, f7.5,1x, es8.2,1x, f7.5,1x, f7.5,1x, es11.5)
!    write(logstring,*) " cycle      NHT      k      time     tacc   nc  kcmfd   psidiff  avgk    stdevk    SHN_E"
10  format     (i6,1x, i10,1x,  f7.5,1x, f8.3,1x, f8.3,1x, i3,1x, f7.5,1x, es8.2,1x, es8.2,1x, f7.5,1x, f7.5,1x, es11.5)
101 format     (i6,1x, i10,1x,  f7.5,1x, f8.3,1x, f8.3,1x, i3,1x, f7.5,1x, es8.2,1x, es8.2,1x, es8.2,1x, f7.5,1x, f7.5,1x, es11.5,1x, es11.5,1x, es11.5,1x, es11.4,1x, es11.4,1x, es11.4,1x, es11.4,1x)
    if (ireft .eq. 1) then
        write(logstring,*) " cycle      NHT      k      time     tacc   nc  kcmfd  err_psi  err_asmp err_p2bf  avgk   stdevk     SHN_E"
    else
        write(logstring,*) " cycle      NHT      k      time     tacc   nc  kcmfd   psidiff  n_balance avgk    stdevk    SHN_E"
    endif
11  format     (i6,1x, i10,1x,  f7.5,1x, f8.3,1x, f8.3,1x, f7.5,1x, f7.5,1x, es11.5)
    call log_(logstring)
    
12  format (100(es11.4, 1x))

    nwt = history
    if (.not.(lact)) nwt=nwt*mltfinact 
    

    nrmj=1._8/hcmx
    nrmp=1._8
    nrmpsi=1._8/vfm
    vcmr=1._8/vcm  
    scmrx=1._8/hcmx
    scmry=1._8/hcmy
    
    
    ke = 0.0_8
    tmcacc = 0.0_8    
    lmwnd = .false.
    call init_statistics
    
    if (lact .eq. .false.)then
!        if (CMFDIACT) bFeedback = .true.
        k_nn = kinit
        k_mc = kinit
        k_cmfd = kinit
        k_cmfdp = k_cmfd
!    else
!        if (CMFDACT) bFeedback = .true.
    endif
    
    
    !call make_initial_particle(nq, ntrn, bat)
    !call initp_from_dist(nq, ntrn, bat)
    if (lact .eq. .false. .and. swfixtype .eq. 0) then
        nq = nwt
        if (lmgfdmfirst .or. lreadFDM)then
            call FS_FDM(nq, nqtot, nwt, ntrn(:))
    		nsource = nqtot
        else
            call initp_from_dist(nq, ntrn)
    		nsource = nwt
        endif
        gseed=strideN(gseed, nqtot)

    else
        
    endif
    
#ifdef _OPENMP
!$  ts(1) = omp_get_wtime()
#else
    call cpu_time(ts(1))
#endif


    do batch = 1, nbatch
        
        ! initialization of local seed
#ifdef MPIF
        if (myrank>0) lseed=strideN(gseed, sum(nqmpi(:,myrank-1)))
#else
        lseed=gseed
#endif
        
        testnnt=0
        mccount = mccount+1
        ! check CMFD condition
        call cmfd_ctrl(lact, batch, loncmfd, lonfdb) 

	    gk = 0.
        GCTphiv=0;GCTtot=0;GCTnuf=0;GCTabsr=0;GCTsct=0;
        phisuu=0;phiau=0;
!        tphisuu=0;tphiau=0;
        psisumacc=0;
        psisumcyc=0;

		maxpinerr =0;
        tphi2x=0._8
        tphi2y=0._8
        tphi2trx=0._8
        tphi2try=0._8
        ttrcwx=0._8
        ttrcwy=0._8
        
#ifdef ASMPRL            
        tnfs=0
        tnin=0
        tnout=0
        tnabsr=0
#endif        
        
        fis1_4=0
        fis1_4acc=0

!$        nhtcnt = 0
!$        gkomp  = 0

        
        if (delkw>0._8) then
            ke = k_cmfd + shiftk(batch, delkw)
            k_nn = 1._8/(1._8/k_nn-1._8/ke)
        endif  
        ntally = 0
        nc_super = 1
        nhtcnt = 0
        tid = 1
820		format (i2, 1x, 32(i4, 1x))		
!		print 820,myrank, ntrnomp(1,1,1:8)%lm(1)
!!$      nqomp = 0
!$omp   parallel private(tid, mi, nnt, lkomp, nqompt, base, g, g2, i, j, tpin, nr, l, r, k, t, jc, ic, gp, volflux, phir, rphiv, qi, qj, tpsi, pinerr)&
!$omp   default(shared)
!$      tid=omp_get_thread_num()+1

! =======================================================================
! ======== Initialize variables =========================================
!$      lkomp = 0
!$      nqompt = 0
!$      jcomp(:,:,:,:,tid)=0
!$      psiomp(:,:,tid)=0;
!$      nfsrcth(:,:,tid)=0;
!!$      nfsrcfth(:,:,tid)=0;
!$omp barrier
! ======== Initialize variables =========================================
! =======================================================================
    
        do while(nq>0)
        nhtcnt(tid) = 0
!$omp do reduction(+:nhtcnt)
            do mi=1, nq
!                write(*,*) i
#ifdef _OPENMP      
                nhtcnt(tid) = nhtcnt(tid)+1
                call migrate(ntrn(mi), k_nn*swsplitting, ke*swsplitting, nnt, lseed, mi, &
                tid, psiomp(:,:,tid), jcomp(:,:,:,:,tid), gkomp(tid), ntrnomp(:,:,tid), nqomp(:,tid))
#else
                call migrate(ntrn(mi,bat(NOW)), k_nn*swsplitting, ke*swsplitting, nnt)
#endif
            end do
!$omp end do
!$omp barrier ! wielandt or super history method treatment
            base=sum(nqomp(OMPCUR,1:tid-1))
            do i=1, nqomp(OMPCUR,tid)
                ntrn(base+i)=ntrnomp(i,OMPCUR,tid)
            enddo
!$omp barrier
!$omp master
            ntally = ntally + nq
            nq = sum(nqomp(OMPCUR,:))
            nqomp(OMPCUR,:)=0
!$omp end master
!$omp barrier
            if (SUPERHST>0)then
                if (nc_super .eq. SUPERHST) then
                    exit
                else
!                    call switch_cycle(bat)
                    nc_super  = nc_super +1
                endif
            endif 
!            if (ke>0._8 .and. SUPERHST .eq. 0) call switch_wielandt(bat)
        end do

!$omp barrier
!        print*, tid, nqomp(OMPNXT, tid), sum(nfsrcth(:,:,tid))
        base=sum(nqomp(OMPNXT,1:tid-1))
        do i=1, nqomp(OMPNXT,tid)
            ntrn(base+i)=ntrnomp(i,OMPNXT,tid)
        enddo
!$omp barrier
!!$omp master
!!	do j=1, nthrd
!!		base=sum(nqomp(OMPNXT,1:j-1))
!!		if (myrank .eq. 0) print *, base
!!		do i=1, nqomp(OMPNXT,j)
!!			ntrn(base+i)=ntrnomp(i,OMPNXT,j)
!!		enddo
!!	enddo
!!$omp end master
!$omp master
        nq = sum(nqomp(OMPNXT,:))
        nqtot=nq
        nqomp(OMPNXT,:)=0

        gk=sum(gkomp)
!$omp end master

        ! collect fission source distribution
        lsupacc(tid)=0;lsupcyc(tid)=0;lslw(tid)=0
    
!$omp barrier
!$omp do    
    do j=1, nfmy;do i=1, nfmx
!        acpsip(i,j)=acpsi(i,j)
        psifp(i,j)=psif(i,j)
        tpsi = sum(psiomp(i,j,:))*nrmpsi*volfr_r(i,j)
        psif(i,j)=tpsi
!        nfsrcf(i,j)=nfsrcf(i,j)+sum(nfsrcfth(i,j,:))
    enddo;enddo
!$omp end do  

#ifdef ASMPRL            
!$omp do    
    do j=1, ncmy;do i=1, ncmx
        nfs(i,j)=nfs(i,j)+sum(tnfs(i,j,:))
        nin(i,j)=nin(i,j)+sum(tnin(i,j,:))
        nout(i,j)=nout(i,j)+sum(tnout(i,j,:))
        nabsr(i,j)=nabsr(i,j)+sum(tnabsr(i,j,:))
        nfsrc(i,j)=nfsrc(i,j)+sum(nfsrcth(i,j,:))
    enddo;enddo
!$omp end do     
#endif
!$omp barrier

!=======================================================================
!$omp barrier
!=======================================================================
#ifndef _NOCMFD
!=======================================================================
!===================== homogenization ==================================

! homogenization in MC simulation
    if( ((CMFDIACT .eq. .true.) .and. (lact .eq. .false.)) .or. ((CMFDACT .eq. .true.) .and. (lact .eq. .true.)) ) then ! .or. (CMFDACT .and. lact ))
    if (lnoacc) then
!$omp do
    ! summation
    do j=1,ncmy;do g=1,ngcmfd;do i=1,ncmx
        GCAphiv(g,i,j)  =0._8
        GCAtot(g,i,j)   =0._8
        GCAnuf(g,i,j)   =0._8
        GCAabsr(g,i,j)  =0._8
        GCAsct(:,g,i,j) =0._8
        hjac(:,g,i,j)   =0._8
        
        trcwx(g,i,j)    =0._8
        trcwy(g,i,j)    =0._8
    enddo;enddo;enddo
!$omp enddo
!$omp do
    do j=0,ncmy;do g=1,ngcmfd;do i=0,ncmx
        phi2x(g,i,j)  =0._8
        phi2y(g,i,j)  =0._8
        phi2trx(g,i,j)=0._8
        phi2try(g,i,j)=0._8
    enddo;enddo;enddo
!$omp enddo
    endif
!$omp do
    ! summation
    do j=1,ncmy;do i=1,ncmx;do g=1,ngcmfd
        GCAphiv(g,i,j)  =GCAphiv(g,i,j) +sum(GCTphiv(g,i,j,:))
        GCAtot(g,i,j)   =GCAtot(g,i,j)  +sum(GCTtot(g,i,j,:))
        GCAnuf(g,i,j)   =GCAnuf(g,i,j)  +sum(GCTnuf(g,i,j,:))
        GCAabsr(g,i,j)  =GCAabsr(g,i,j) +sum(GCTabsr(g,i,j,:))
        do gp=1,ngcmfd
            GCAsct(gp,g,i,j) =GCAsct(gp,g,i,j)+sum(GCTsct(gp,g,i,j,:))
        enddo
        ! GCMFD
        trcwx(g,i,j)=trcwx(g,i,j)+sum(ttrcwx(g,i,j,:))*vcmr
        trcwy(g,i,j)=trcwy(g,i,j)+sum(ttrcwy(g,i,j,:))*vcmr
    enddo;enddo;enddo
!$omp enddo


!$omp do
    do j=0,ncmy;do i=0,ncmx;do g=1,ngcmfd
        phi2x(g,i,j)  =phi2x(g,i,j)  +sum(tphi2x(g,i,j,:))*scmry
        phi2y(g,i,j)  =phi2y(g,i,j)  +sum(tphi2y(g,i,j,:))*scmrx
        phi2trx(g,i,j)=phi2trx(g,i,j)+sum(tphi2trx(g,i,j,:))*scmrx
        phi2try(g,i,j)=phi2try(g,i,j)+sum(tphi2try(g,i,j,:))*scmry
    enddo;enddo;enddo
!$omp enddo



!! for cycle CMFD, hj and phivac should be set to zero.
!$omp do
    do j=1, ncmy;do l=1, nthrd;do i=1, ncmx;do g=1, ngcmfd;do k=1, 4
        hjac(k,g,i,j)=hjac(k,g,i,j)+jcomp(k,g,i,j,l)*nrmj
    enddo;enddo;enddo;enddo;enddo
!$omp enddo

    endif
#endif
!$omp end parallel

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Communication throught MPI
#ifdef MPIF	
	call reduce(gk, gkrcv, mpi_mc_comm, .true.)
#ifdef MPIDEBUG
!	if(myrank .eq. 0 )	print *, myrank, gk, gkrcv
#endif
	gk=gkrcv

	call reduce(psif, psifrcv, nfmx, nfmy, mpi_mc_comm, .false.)
	call reduce(GCAtot, totrcv, ngcmfd, ncmx, ncmy, mpi_mc_comm, .false.)
	call reduce(GCAnuf, nufrcv, ngcmfd, ncmx, ncmy, mpi_mc_comm, .false.)
	call reduce(GCAphiv, phivrcv, ngcmfd, ncmx, ncmy, mpi_mc_comm, .false.)
	call reduce(GCAsct, sctrcv,ngcmfd, ngcmfd, ncmx, ncmy, mpi_mc_comm, .false.)
	call reduce(hjac, jcrcv, 4, ngcmfd, ncmx+2, ncmy+2, mpi_mc_comm, .false.)
!	call testnn(ntrn(:, bat(NXT)), nq(bat(NXT)), nqtot)
	call nbalance(ntrn(:), nq, nqtot)

	! for GCMFD
	call reduce(trcwx, trcwxrcv, ngcmfd, ncmx, ncmy, mpi_mc_comm, .false.)
	call reduce(trcwy, trcwyrcv, ngcmfd, ncmx, ncmy, mpi_mc_comm, .false.)

	call reduce(phi2x, phi2xrcv, ngcmfd, ncmx+1, ncmy+1, mpi_mc_comm, .false.)
	call reduce(phi2y, phi2yrcv, ngcmfd, ncmx+1, ncmy+1, mpi_mc_comm, .false.)

	call reduce(phi2trx, phi2trxrcv, ngcmfd, ncmx+1, ncmy+1, mpi_mc_comm, .false.)
	call reduce(phi2try, phi2tryrcv, ngcmfd, ncmx+1, ncmy+1, mpi_mc_comm, .false.)

	trcwx=trcwxrcv
	trcwy=trcwyrcv
	phi2x=phi2xrcv
	phi2y=phi2yrcv
	phi2trx=phi2trxrcv
	phi2try=phi2tryrcv
#ifdef MPIDEBUG
!	print *, "Source Neutron in Queue", myrank, nq(bat(nxt)), nqtot
#endif

#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


    
!    psia=0._8;

!$omp parallel private(tid, mi, nnt, lseed, lkomp, nqompt, base, g, g2, i, j, tpin, nr, l, r, k, t, jc, ic, gp, volflux, phir, rphiv, qi, qj, tpsi, ja, ia, sty, stx)&
!$omp default (shared)
!$      tid=omp_get_thread_num()+1
! MC tallies ==================================================================================
!$omp do    
    do j=1, nfmy;ja=f2cy(j);do i=1, nfmx;ia=f2cx(i)
        !tpsi=ppsif(i,j)
        tpsi=ppsif(i,j)
        stpsi(i,j)=stpsi(i,j)+tpsi
        stpsisq(i,j)=stpsisq(i,j)+tpsi**2
        psisumcyc(tid)=psisumcyc(tid)+tpsi
        psisumacc(tid)=psisumacc(tid)+stpsi(i,j)
        
        ! check tilt
        qi=qtmapx(i); qj=qtmapy(j)
        fis1_4(qi,qj,tid)=fis1_4(qi,qj,tid)+psif(i,j)
        fis1_4acc(qi,qj,tid)=fis1_4acc(qi,qj,tid)+stpsi(i,j)        
    enddo;enddo
!$omp end do
!$omp do    
    do j=1, nfmy
        !tpsi=ppsif(i,j)
        tpsi=sum(ppsif(:,j))
        stpsiz(j)=stpsiz(j)+tpsi
        stpsizsq(j)=stpsizsq(j)+tpsi**2
        ! check tilt
    enddo
!$omp end do!$omp barrier
!    if (not(hybridcmfd )) then
!$omp do  ! loop for assembly power
    do ja=1, ncmy; sty=(ja-1)*rcy; do ia=1, ncmx; stx=(ia-1)*rcx
        l=(ja-1)*ncmx+ia
        tpsi=sum(ppsif(stx+1:stx+rcx, sty+1:sty+rcy))
        stpsia(l,1)=stpsia(l,1)+tpsi
        stpsiasq(l,1)=stpsiasq(l,1)+tpsi**2
    enddo; enddo
!$omp enddo
!    endif
!$omp barrier
!$omp single 
    rnormpsiacc=1._8/sum(psisumacc)
    rnormpsicyc=1._8/sum(psisumcyc)
!$omp end single
!$omp barrier
    if (ireft .eq. 1) then
!$omp do    
        do j=1, nfmy;do i=1, nfmx
            lsupacc(tid)=lsupacc(tid)+(stpsi(i,j)*rnormpsiacc-refpsi(i,j))**2
            !pinerr=abs(stpsi(i,j)*rnormpsiacc-refpsi(i,j))/refpsi(i,j)
            pinerr=abs(ppsif(i,j)*rnormpsicyc-refpsi(i,j))/refpsi(i,j)
            if (refpsi(i,j)*npowpinref>0.5_8 .and. pinerr > maxpinerr(tid)) maxpinerr(tid)=pinerr
            lsupcyc(tid)=lsupcyc(tid)+(ppsif(i,j)*rnormpsicyc-refpsi(i,j))**2
        enddo;enddo
!$omp end do    
    else
!$omp do    
        do j=1, nfmy;do i=1, nfmx
            lsupcyc(tid)=lsupcyc(tid)+(psifp(i,j)-psif(i,j))**2
        enddo;enddo
!$omp end do    
    endif
!$omp barrier
!$omp single ! check 
    if (ireft .eq. 1) then
        perracc = sqrt(sum(lsupacc(:))*npowpinref) !/pinpowsum
        perrcyc = sqrt(sum(lsupcyc(:))*npowpinref) !/pinpowsum
    else
        fpsidiff = sqrt(sum(lsupcyc(:))*npowpinref)/sum(psisumcyc)
    endif
    do i=2, nthrd
        if (maxpinerr(i)>maxpinerr(1)) maxpinerr(1)=maxpinerr(i)
    enddo
!$omp end single
! =============================================================================================

#ifndef _NOCMFD
#ifdef MPIF
	if (myrank .eq. 0) then
#endif
!=======================================================================
!===================== homogenization ==================================

! PREPARE TALLY DISCARD IN THE LAST CMFD.
! reset reaction rates.
!$omp barrier
!$omp single
    if (.not.(RSTDHAT) .and. CMFDDISCARD) then
        if (mccount .eq. ndscrd) then
            rjds    =pjc   
            rphids  =pphiv
            rtotds  =ptot 
            rNuFds  =pnuf 
            rmsctds =psct 
        elseif (mccount .eq. (inactive+active))then
            pjc  =pjc  -rjds   
            pphiv=pphiv-rphids 
            ptot =ptot -rtotds 
            pnuf =pnuf -rNuFds 
            psct =psct -rmsctds
            print *, "Tallies are discarded for the last CMFD calculation.", ndscrd
        endif
    endif
!$omp end single
!$omp barrier

! homogenization in MC simulation
    if( ((CMFDIACT .eq. .true.) .and. (lact .eq. .false.)) .or. ((CMFDACT .eq. .true.) .and. (lact .eq. .true.)) ) then ! .or. (CMFDACT .and. lact ))

!$omp barrier
!$omp do
    ! generate group constant
    do j=1, ncmy;do i=1, ncmx; do g=1, ngcmfd
        if (pphiv(g,i,j)>0) then
            rphiv=1._8/pphiv(g,i,j)
            htot(g,i,j)     =ptot(g,i,j)*rphiv
            hnuf(g,i,j)     =pnuf(g,i,j)*rphiv
            hphi(g,i,j)     =pphiv(g,i,j)*vcmr
            hdif(g,i,j)     =1._8/3/hTot(g,i,j)
            do gp=1, ngcmfd
                hmsct(gp,g,i,j)  =psct(gp,g,i,j)*rphiv
            enddo            
            hRmv(g,i,j)     =htot(g,i,j)-hmsct(g,g,i,j)
        else
            htot(g,i,j)     =0._8
            hnuf(g,i,j)     =0._8
            hphi(g,i,j)     =0._8
            hdif(g,i,j)     =0._8
            hmsct(:,g,i,j)  =0._8
            hRmv(g,i,j)     =0._8
        endif
    enddo;enddo;enddo    
!$omp enddo
!$omp barrier
!$omp do
    do j=1, ncmy;do i=1, ncmx;do g=1, ngcmfd;do k=1, 4
        hj(k,g,i,j)=pjc(k,g,i,j)-pjc(LTOPST(k),g,i+LTNEIGH(1,k),j+LTNEIGH(2,k))
    enddo;enddo;enddo;enddo
!$omp enddo
    endif
#ifdef MPIF
	endif 
#endif
#endif
!$omp end parallel
    
    ! pin2box factor
    powasm=0.
    p2b=0.
    SHNp2bc=0.
    do j=1, ncmy
        aj=(j-1)*rcy+1
        do i=1, ncmx
            ai=(i-1)*rcx+1
            powcm(i,j)=sum(ppsif(ai:ai+rcx-1,aj:aj+rcy-1))
            if (powcm(i,j)>0) then
                do l=1, rcy
                    fj=aj+l-1
                    do k=1, rcx
                        fi=ai+k-1
                        if (ppsif(fi, fj)>0)then
                            p2b(k,l,i,j) = ppsif(fi,fj)/powcm(i,j)
                            SHNp2bc(i,j) = SHNp2bc(i,j)-p2b(k,l,i,j)*log(p2b(k,l,i,j))
                        endif
                    enddo
                enddo
                SHNp2b = SHNp2bc(i,j)*powcm(i,j)
            endif
        enddo
    enddo
    SHNp2b = SHNp2b/sum(powcm)
    

#ifndef MPIF
    ! FSD write
    ! check for cycle by cycle fission source distribution
    if (Woffset>0 .and. mod(batch-1,Woffset) .eq. 0) then
        if (and(WFSDinact, not(lact)) .or. and(WFSDact,lact)) then
            write (sfile, '("psi_step",L1,I<MAX(0,0)>,".bin")'), lact,batch
            open(fpsi, file =sfile, form='binary', status='unknown')
            if  (WFSDacc .eq. 1) then
                write(fpsi) stpsi
            else 
                write(fpsi) psif
            endif
            close(fpsi)
        endif        
    endif

        ! tilt check
        tilt_c=0
        tilt_a=0
        do j=1, 2
            do i=1,2
                qtp(i,j)=sum(fis1_4(i,j,:))
                qtpa(i,j)=sum(fis1_4acc(i,j,:))
                ! half of interface
                qtp(i,j)=qtp(i,j)+sum(fis1_4(0,j,:))/2+sum(fis1_4(i,0,:))/2
                qtpa(i,j)=qtpa(i,j)+sum(fis1_4acc(0,j,:))/2+sum(fis1_4acc(i,0,:))/2
            enddo
        enddo
        fss=sum(qtp)
        fssa=sum(qtpa)
        qtp=qtp/sum(qtp)
        qtpa=qtpa/sum(qtpa)
        tiltx = (qtp(2,1)+qtp(2,2)-qtp(1,1)-qtp(1,2))/2
        tilty = (qtp(1,2)+qtp(2,2)-qtp(1,1)-qtp(2,1))/2
        
        tilt_c=0;tilt_a=0
!        do j=1, 2
!            do i=1, 2
!                tilt_c=tilt_c-log(qtp(i,j))
!                tilt_a=tilt_a-log(qtpa(i,j))
!            enddo
!        enddo
!        tilt_c = tilt_c/log(2._8)-8
!        tilt_a = tilt_a/log(2._8)-8

        do j=1, 2
            do i=1, 2
                tilt_c=tilt_c+(1._8/4-qtp(i,j))**2
                tilt_a=tilt_a+(1._8/4-qtpa(i,j))**2
            enddo
        enddo
        tilt_c = sqrt(tilt_c)*sqrt(4._8)
        tilt_a = sqrt(tilt_a)*sqrt(4._8)
        
    !    print 12, tilt_c, tilt_a, tiltx, tilty
    
        if (ireft .eq. 1) then
            ! AUG 18 2010 LMJ
            ! Use reference of coarse mesh power distribution
            ! so only effective when the coarse mesh pitch = assembly pitch
            npxm=npx-1;npxp=npx+1
            npym=npy-1;npyp=npy+1
            ! generate assembly-wise power
            upper=0._8
            do j=1, nay
                aj=(j-1)*npy+1
                do i=1, nax
                    ai=(i-1)*npx+1
                    powasm(i,j)=sum(stpsi(ai:ai+npxm,aj:aj+npym))
                    do l=0, npym
                        fj=aj+l
                        do k=0, npxm
                            fi=ai+k
                            if (stpsi(fi, fj)>0)then
                                upper = upper+ (stpsi(fi,fj)/powasm(i,j)-refp2bf(fi,fj))**2
                            endif
                        enddo
                    enddo
                enddo
            enddo
            powasm=powasm/sum(powasm)
            
            rnormpsiacc=rnormpsiacc*npowpinref/npowasmref
            upper =0
            do j=1, nay
                aj=(j-1)*npy+1
                do i=1, nax
                    ai=(i-1)*npx+1
                    do l=0, npym
                        fj=aj+l
                        do k=0, npxm
                            fi=ai+k
                            if (stpsi(fi, fj)>0)then
                                pin2bf(fi,fj)=stpsi(fi,fj)*rnormpsiacc/powasm(i,j)
                                upper = upper+(pin2bf(fi,fj)-refp2bf(fi,fj))**2
                            endif
                        enddo
                    enddo
                enddo
            enddo           
            err_p2bf=sqrt(upper*npowpinref)/sump2bf
!            print*, sum(pin2bf)
            
            ! check the error of assembly power
            
!            err_asmp=0
!            do j=1, nay
!                do i=1, nax
!                    err_asmp = err_asmp+(trefcmfd(i,j)-powasm(i,j))**2
!                enddo
!            enddo
!            err_asmp=sqrt(err_asmp*npowasmref)
            
        endif
#endif
        
        !ntallywt = history - nsource + ntally ! weight sum of all neutrons
        call estimateK(nwt, k_mc)
!====================================================================
#ifndef MPIF
! neutron balance check
        errblc = 0
        if (k_cmfd .eq. 0) k_cmfd = k_mc
        if (loncmfd) then
            upper =0; lower=0
            do j=1, ncmy;do i=1, ncmx;do g=1, 2
            
                leak = sum(hj(1:2,g,i,j))*ncmy+sum(hj(3:4,g,i,j))*ncmx
                nrmv = hRmv(g,i,j)*hphi(g,i,j)*vcm

                nfis = 0
                if (g .eq. 1) then
                    nfis = (hnuf(1,i,j)*hphi(1,i,j)+hnuf(2,i,j)*hphi(2,i,j))*vcm/k_cmfd
                endif
                
                nsct = hmSct(g,1,i,j)*hphi(1,i,j)*vcm + hmSct(g,2,i,j)*hphi(2,i,j)*vcm
                nsct = nsct - hmSct(g,g,i,j)*hphi(g,i,j)*vcm 

                lhs = leak + nrmv
                rhs = nfis + nsct
                diff = lhs - rhs
                lower = lower + rhs
                upper = upper + diff**2
                if (g .eq. 1) then
    !                print*, diff
!                    print *, (nsct+diff)/hphi(g,i,j)/vcm, hmSct(1,2,i,j)
!                    hmSct(1,2,i,j)=(nsct+diff)/hphi(g,i,j)/vcm
                endif
                
            enddo;enddo;enddo
            errblc=sqrt(upper)/lower
        endif
#endif
!====================================================================


!#ifdef MPIF
!	if (myrank .eq. 0) then
!#endif
        shnEtrp = ShannonEntropy(nfmx,nfmy,ppsif)
        shnEacc = ShannonEntropy(nfmx,nfmy,stpsi)
        

        ! k_cmfd = k_mc
        NMVWND =1

#ifndef _NOCMFD        
!        call cmfd(k_cmfd, err, iter, FALSE)
#ifdef MPIF
	if (myrank .eq. 0) then
#endif
        iter = 0
        k_cmfd = k_cmfdp
        if (loncmfd) then
        
            if (WBINCMFD) then
                call write_binarydata4CMFD
             endif
            
		!	print *, hrmv(:,:,10)
		!	print *, hmsct(:,:,:,10)
		!	print *, hphi(:,:,10)
		    if (k_cmfd>0 .and. k_cmfd<2.0) then
		    else
		        k_cmfd = k_mc
		    endif
		    if (rungcmfd) then
                call lcmfd2gacc(k_cmfd, err, iter, err_cmfd, lonfdb, powc)
		    else
                call cmfd2gacc(k_cmfd, err, iter, err_cmfd, err_cmfd_max, lonfdb, powc)
		    endif
            
		    if (k_cmfd>0 .and. k_cmfd<2.0) then
		    else
		        k_cmfd = 0._8
		        fdbsafe = .false.
		    endif
            if (fdbsafe) then
                shnECMFD = ShannonEntropy(ncmx, ncmy, adjp)
                if (PPRECONST) then
                    call pinpower_reconstruction(p2b, powc,recpow)
                    shnREC = ShannonEntropy(nfmx, nfmy, recpow)
                    shnECMFD = shnREC
                    rmserr_rec = rmserr(recpow)
                endif
            endif
            
            if (lfirstcmfd) then ! reset accumulated tallies
!                print *, "RESET"
                lfirstcmfd = .false. 
!                do g=1, ng;do j=1, nfmy;do i=1, nfmx
!!                    jfAC(:,i,j,g) = 0.
!!                    phivac(:,i,j,g) = 0.
!                enddo;enddo;enddo
            endif
            k_cmfdp = k_cmfd
			if (nignracc >= batch .and. lact .eq. .false.) then
!$omp parallel private(i,j,g)
!$omp do
				! tally initialization
				do j=1,ncmy;do i=1,ncmx;do g=1,ngcmfd
					GCAphiv(g,i,j)  =0._8
					GCAtot(g,i,j)   =0._8
					GCAnuf(g,i,j)   =0._8
					GCAabsr(g,i,j)  =0._8
					GCAsct(:,g,i,j) =0._8
					! GCMFD
					trcwx(g,i,j) = 0._8
					trcwy(g,i,j) = 0._8
				enddo;enddo;enddo
!$omp enddo

!$omp do
				do j=0,ncmy;do i=0,ncmx;do g=1,ngcmfd
					phi2x(g,i,j)   = 0._8
					phi2y(g,i,j)   = 0._8
					phi2trx(g,i,j) = 0._8
					phi2try(g,i,j) = 0._8
				enddo;enddo;enddo
!$omp enddo
!$omp end parallel
			endif
		else 
            k_cmfd =-0._8
        endif
        
!        call solFMC
        call add_statistics(k_mc, k_cmfd)
#ifdef MPIF
	endif
#endif

#endif
#ifdef MPIF
	if (myrank .eq. 0) then
#endif
#ifdef _OPENMP
!$  ts(2) = omp_get_wtime()
#else
    call cpu_time(ts(2))
#endif

    gseed=strideN(gseed, nsource)

    tmc = ts(2) - ts(1)
    tmcacc = tmcacc + tmc
    ts(1)=ts(2)
#ifdef _NOCMFD
    write(logstring,11) batch, nsource, k_mc, tmc, tmcacc, avgk_mc, stdk_mc, shnEtrp
#else
    if (loncmfd) then
!        write(logstring,10) batch, nsource, k_mc, tmc, tmcacc, iter, k_cmfd, fpsidiff, err_cmfd, avgk_mc, stdk_mc, shnEtrp
        if (ireft .eq. 1) then 
            call find_errdist(s1, s2)
            if (ppreconst) then
                write(logstring,101) batch, nsource, k_mc, tmc, tmcacc, iter, k_cmfd, fpsidiff, err_asmp, err_p2bf, avgk_mc, stdk_mc, shnEtrp, shnEAcc, shnECMFD, perrcyc,   rmserr_rec, err_cmfd_max,  s2 !maxpinerr(1)
            else
                write(logstring,101) batch, nsource, k_mc, tmc, tmcacc, iter, k_cmfd, fpsidiff, err_asmp, err_p2bf, avgk_mc, stdk_mc, shnEtrp, shnEAcc, shnECMFD, perrcyc,   perracc, err_cmfd, s2 !maxpinerr(1)
            endif
!            write(logstring,101) batch, nsource, k_mc, tmc, tmcacc, iter, k_cmfd, fpsidiff, err_asmp, err_p2bf, avgk_mc, stdk_mc, shnEtrp, SHNp2b, shnECMFD, perrcyc,   perracc, err_cmfd, maxpinerr(1)
        else
            write(logstring,101) batch, nsource, k_mc, tmc, tmcacc, iter, k_cmfd, fpsidiff, err_asmp, err_p2bf, avgk_mc, stdk_mc, shnEtrp, shnEAcc, shnECMFD, tiltx, tilty, gtally(1), 0
!            write(logstring,101) batch, nsource, k_mc, tmc, tmcacc, iter, k_cmfd, fpsidiff, err_asmp, err_p2bf, avgk_mc, stdk_mc, shnEtrp, SHNp2b, shnECMFD, tiltx, tilty, gtally(1), 0
        endif
    else
        write(logstring,10) batch, nsource, k_mc, tmc, tmcacc, iter, k_cmfd, perrcyc, perracc, avgk_mc, stdk_mc, shnEtrp
    endif
!    print*, fpsidiff
#endif
    call log_(logstring) 
#ifdef MPIF  ! the end of subroutines for master processor
	endif
	call bcast(fdbsafe, mpi_mc_comm)
#endif
        if (lastcycle_nshrink()) then
            nq = history
            nwt = history
            nsource = history
        endif
        
        !** process for next step **!
        if (lonfdb .and. fdbsafe) then
            if (lresampling .and. rcx .eq. 1 .and. rcy .eq. 1) then
                call FS_CMFD(nq, nqtot, nwt, ntrn(:))
            else
                call feedback(nq, nwt, ntrn(:))
            endif
        else
            call init_weight(nq, nwt, ntrn(:))
        end if
		if (.not.(lastcycle_nshrink())) nsource = nqtot 

!        k_nn = k_mc
        if ((swfeedback .eq. .false.) .and. (swfixtype .eq. 0)) then
!            k_nn = k_mc*nqn/history
            k_nn = k_mc
!             k_nn = k_mc
        else 
!            k_nn = k_mc
            k_nn = k_mc
        endif
        
!        call switch_cycle(bat)
        
        if (WRTREF) then
            open(2001, file='REF.bin', form='binary', status='unknown')
            write(2001) phivac
            close(2001)
        endif 
        totbat = totbat+1
#ifdef _THFDB
        call fdb_th  
#endif

    end do
    contains
    
subroutine find_errdist(s1, s2)
    use ifport
    include "arrays.FI"
    real(8), pointer, dimension(:,:) :: sol
    real(8), intent(inout) :: s1, s2
    real(8) :: sumsol, norm
    integer :: l, aim1, aim2
    real(8) :: cmp1, cmp2
    integer :: ret
    real(8), dimension(nfmx) :: errline
    integer, save :: first = .true.
    
    if (first) then
        ret = MAKEDIRQQ("errdist")
        first = .false.
    endif
    if (PPRECONST) then
        ! final solution is recpow
        sol=>recpow
    else
        ! final solution is psiacc
        sol=>stpsi
    endif
    
    sumsol=sum(sol)
    s1=0;s2=0
    if (sumsol>0) then
        if (WERRDIST) then
            write (sfile, '("./errdist/err",I<MAX(0,0)>,".txt")'), batch
            open (123, file = sfile, status='unknown')
        endif
        norm = 1._8/sumsol
        l=1
        do j=1, nfmy
            errline=0
            do i=1, nfmx
                if (refpsi(i,j) .eq. 0) cycle
                errline(i) =(sol(i,j)*norm-refpsi(i,j))/refpsi(i,j)
                errdist(l)=abs(errline(i))
                l=l+1
            enddo
                if (WERRDIST) write(123, '(2000(es12.5,1x))'), errline
        enddo
        if (WERRDIST) close(123)
        aim1 = npowpinref*0.68268949_8
        aim2 = npowpinref*0.954499736_8
!        aim1 = npowpinref-1
!        aim2 = npowpinref-2
        
!        s1=quicksearch(errdist,1,npowpinref, aim1)
        s2=quicksearch(errdist,1,npowpinref, aim2)
    endif
endsubroutine    
    
function lastcycle_nshrink 
    logical lastcycle_nshrink 
    lastcycle_nshrink = .false.
    if (batch .eq. nbatch) then
        if (.not.(lact)) then
            if (mltfinact>1.0) then
                lastcycle_nshrink = .true.
            endif
        endif
    endif
end function

subroutine pinpower_reconstruction(p2bc, asmp, powout)
    include "pspec.FI"
    integer, save :: called = 0 
    real(8), save, pointer, dimension(:,:,:,:)  :: p2bacc
    real(8), dimension(:,:,:,:), intent(in) :: p2bc
    real(8), dimension(:,:), intent(in) :: asmp
    real(8), dimension(:,:), intent(out) ::powout
    integer :: i, j, l, k, fi, fj, is, ie, js, je
    real(8) :: rsum, rnorm
    
    called = called+1
    if (called .eq. 1) then
        allocate(p2bacc(rcx, rcy, ncmx, ncmy))
        p2bacc = 0
    endif
    if (called > nskip_p2b) then
        do j=1, ncmy; do i=1, ncmx
            p2bacc(:,:,i,j)=p2bacc(:,:,i,j)+p2bc(:,:,i,j)
        enddo;enddo
        fi=0;fj=0;
        do j=1, ncmy; js=(j-1)*rcy+1;je=j*rcy;
        do i=1, ncmx; is=(i-1)*rcx+1;ie=i*rcx;
            rsum = sum(p2bacc(:,:,i,j))
            if (rsum>0) then
                rnorm = asmp(i,j)/rsum
                fj=js
                do k=1, rcy;fi=is
                do l=1, rcx
                    powout(fi,fj) = p2bacc(l,k,i,j)*rnorm
                    fi=fi+1
                enddo;fj=fj+1
                enddo
            endif
        enddo;enddo
    endif
end subroutine

function rmserr(powin)
    real(8), dimension(:,:) :: powin
    integer :: i, j
    real(8) :: rmserr, norm
    norm = sum(powin)
    if (ireft .eq. 1) then
        do j=1, nfmy;do i=1, nfmx
            rmserr=rmserr+(powin(i,j)*norm-refpsi(i,j))**2
        enddo;enddo
    endif
    rmserr = sqrt(rmserr*npowpinref)
end function

subroutine write_binarydata4CMFD
    use ifport
	integer, parameter :: fbincmfd = 1020
	integer :: ret
    integer, save :: first = .true.
    if (first) then
        ret = MAKEDIRQQ("bincmfd")
        first = .false.
    endif

    write (sfile, '("./bincmfd/jc",I<MAX(0,0)>,".bin")'), mccount
    open(fbincmfd, file=sfile, form='binary', status='unknown')
    write(fbincmfd) pjc
    close(fbincmfd)

    write (sfile, '("./bincmfd/phi",I<MAX(0,0)>,".bin")'), mccount
    open(fbincmfd, file=sfile, form='binary', status='unknown')
    write(fbincmfd) pphiv
    close(fbincmfd)

    write (sfile, '("./bincmfd/tot",I<MAX(0,0)>,".bin")'), mccount
    open(fbincmfd, file=sfile, form='binary', status='unknown')
    write(fbincmfd) ptot
    close(fbincmfd)

    write (sfile, '("./bincmfd/nuf",I<MAX(0,0)>,".bin")'), mccount
    open(fbincmfd, file=sfile, form='binary', status='unknown')
    write(fbincmfd) pnuf
    close(fbincmfd)

    write (sfile, '("./bincmfd/sct",I<MAX(0,0)>,".bin")'), mccount
    open(fbincmfd, file=sfile, form='binary', status='unknown')
    write(fbincmfd) psct
    close(fbincmfd)
    
    if (rungcmfd) then
        write (sfile, '("./bincmfd/trcwx",I<MAX(0,0)>,".bin")'), mccount
        open(fbincmfd, file=sfile, form='binary', status='unknown')
        write(fbincmfd) trcwx
        close(fbincmfd)

        write (sfile, '("./bincmfd/trcwy",I<MAX(0,0)>,".bin")'), mccount
        open(fbincmfd, file=sfile, form='binary', status='unknown')
        write(fbincmfd) trcwy
        close(fbincmfd)

        write (sfile, '("./bincmfd/phi2x",I<MAX(0,0)>,".bin")'), mccount
        open(fbincmfd, file=sfile, form='binary', status='unknown')
        write(fbincmfd) phi2x
        close(fbincmfd)

        write (sfile, '("./bincmfd/phi2y",I<MAX(0,0)>,".bin")'), mccount
        open(fbincmfd, file=sfile, form='binary', status='unknown')
        write(fbincmfd) phi2y
        close(fbincmfd)

        write (sfile, '("./bincmfd/phi2trx",I<MAX(0,0)>,".bin")'), mccount
        open(fbincmfd, file=sfile, form='binary', status='unknown')
        write(fbincmfd) phi2trx
        close(fbincmfd)

        write (sfile, '("./bincmfd/phi2try",I<MAX(0,0)>,".bin")'), mccount
        open(fbincmfd, file=sfile, form='binary', status='unknown')
        write(fbincmfd) phi2try
        close(fbincmfd)
    endif
    

end subroutine


end subroutine ! The end of the function MCrun

!subroutine collectntrn
!    use define
!    implicit none
!    include "pspec.FI"
!    include "mcbase.FI"
!    include "mtbuf.FI"
!    integer :: tid
!    integer ::i
!    nqidx = 0
!    do i=1, nthrd-1
!       nqidx(i+1) = nqidx(i)+nqomp(bat(NXT),i)   
!    end do
!    do tid =1, nthrd
!        do i=1, nqomp(bat(NXT),tid)   
!            ntrn(nqidx(tid)+i) = ntrnomp(i, bat(NXT),tid)
!        end do
!    enddo
!!    nq(bat(NXT)) = sum(nqomp(bat(NXT),:))
!end subroutine

subroutine MC(initk)
    use allocs
    use define
    use util_mc
    use RNG
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    include "files.FI"
!$  include "omp_lib.h"    
   
    integer :: nbatch, ichar1st, i
    real(8), intent(in) ::initk
    real(8) :: tchk(2), ct, dummy
    character*200 strseed
    integer :: iseed

    call init_util_mc
    call initRN()
    
    call getarg(1, strseed)
    ichar1st=ichar(strseed(1:1))
    read(strseed, '(i)'), iseed
    call setSEED(getSEED() + (iseed+1)*2**16)
    gseed = getSEED() 

    logstring = "==================================================================================="
    call log_(logstring)
    
!#ifdef RNGMCNP
!    logstring = "RNG : MCNP"
!    call log_(logstring) 
!#endif    
!#ifdef RNGMCCARD
!    logstring = "RNG : MCCARD"
!    call log_(logstring) 
!#endif
#ifdef _IMPCAPT
        logstring = "Implicit capture on"
        call log_(logstring) 
#endif    
#ifdef _EXPTR
        logstring = "Exponential transform on"
        call log_(logstring) 
#endif    

#ifdef _LANISO
    logstring = "Anisotropic scattering with Larsen's Xsection"
    call log_(logstring)
    call init_anisoSct
#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! inactive cycle 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    call prp_inactive
    call MCrun(inactive, initk, .false.)
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! active cycle 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    call prp_active
    call MCrun(active, initk, .true.)

#ifdef MPIF
    if (myrank .eq. 0) &
#endif
    call write_result
end subroutine


subroutine prp_inactive
    use define
    implicit none
    include "pspec.FI"
    include "mcbase.FI"
    include "arrays.FI"
    include "accm.FI"
    include "statistics.FI"
!    character*200 logstring
    
!    bat(NOW) = NOW
!    bat(CUR) = CUR
!    bat(NXT) = NXT

    jfAC     = 0.
    psifAC   = 0.
    phivac   = 0.
    stphi    = 0.
    stphisq  = 0.
    phisAC   = 0.   
    phisaAC  = 0.  
!    trk_estT = 0.

    swfeedback = .false. 
    if (CMFDIACT) swfeedback = .true.
    
#ifdef CMFDTEST
    swfixtype = 0
    swfmshape = 0
    swsplitting = 1
#else
    swfixtype = FXDSRC
    swfmshape = FMSHAPE
    swsplitting = SRCSPLT
#endif

    write(logstring,*) ""
    call log_(logstring)

    logstring = "Inactive cycle"
    call log_(logstring)
    if (swfixtype > 0) then
        write(logstring, '(a,i01)'), "Fixed Source MC in Inactive Cycle  type:",swfixtype 
        call log_(logstring)    
    endif
    
    if (swsplitting .eq. 2) then
        logstring = "Source Spliting"
        call log_(logstring)    
    endif
    if (swfmshape) then
        logstring = "Linear Approximation of flux in a fine mesh"
        call log_(logstring)    
    endif
    if (swfeedback) then
        logstring = "CMFD acceleration on"
        call log_(logstring) 
    else
        logstring = "CMFD acceleration off"
        call log_(logstring) 
    endif
end subroutine


subroutine prp_active
    use lcmfd2g, only: resetpsi_lcmfd2g
    use cmfd2g, only: resetpsi_cmfd2g
    
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "accm.FI"
    include "statistics.FI"
    include "hmgz.FI"
    include "GroupConstant.FI"
!    character*200 logstring
    swfeedback = .false. 
    if (CMFDACT) swfeedback = .true.
#ifndef CMFDTEST
    swfixtype = 0
    swfmshape = 0
    swsplitting = 1
#else
    swfixtype = FXDSRC
    swfmshape = FMSHAPE
    swsplitting = SRCSPLT
#endif

#ifdef ASMPRL  
    nfs=0          
    nin=0
    nout=0
    nabsr=0
#endif

    stphi    = 0.
    stphisq  = 0.
    stpsiz    = 0.
    stpsizsq  = 0.
    stpsi    = 0.
    stpsisq  = 0.
    stpsia    = 0.
    stpsiasq  = 0.
    stpp    = 0.
    stppsq  = 0.
    psifAC   = 0.
!    if (RSTDHAT) then
!        logstring = "Reset Dhat."
!        call log_(logstring)
!!        jfAC = 0.
!!        phivac = 0.
!        hjac=0
!        GCAphiv=0
!        GCAtot=0
!        GCAnuf=0
!        GCAabsr=0
!        GCAsct=0
!
!!        trk_estT = 0.
!    endif
    
    if (rungcmfd) then
        call resetpsi_lcmfd2g()
    else
        call resetpsi_cmfd2g()
    endif

    write(logstring,*) ""
    call log_(logstring)
    logstring = "Active cycle"
    call log_(logstring)
    if (swfixtype > 0) then
        write(logstring, '(a,i01)'), "Fixed Source MC in Inactive Cycle  type:",swfixtype 
        call log_(logstring)    
    endif
    
    if (swsplitting .eq. 2) then
        logstring = "Source Spliting"
        call log_(logstring)    
    endif
    if (swfmshape) then
        logstring = "Linear Approximation of flux in a fine mesh"
        call log_(logstring)    
    endif
    if (swfeedback) then
        logstring = "CMFD acceleration on"
        call log_(logstring) 
    else
        logstring = "CMFD acceleration off"
        call log_(logstring) 
    endif
end subroutine




